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Abstract 

In  this  paper  we  consider  an  application  of  Sobolev-orthogonal  functions  and  radial  basis 
function  to  the  numerical  solution  of  partial  differential  equations.  We  develop  the  funda- 
mentals of  a spectral  method,  present  examples  via  reaction-diffusion  partial  differential 
equations  and  discuss  briefly  some  links  with  theory  of  wavelets. 

1 Introduction 

Radial  basis  functions  are  a well-known  and  useful  tool  for  functional  approximation  in 
one  or  more  dimensions.  The  general  form  of  approximations  is  always  a linear  combin- 
ation (finite  or  infinite)  number  of  shifts  of  a single  function,  the  radial  basis  function. 
In  more  than  one  dimension,  this  function  is  made  rotationally  invariant  by  composing 
a univariate  function,  usually  called  (j>,  with  the  Euclidean  norm.  In  one  dimension  such 
approximation  usually  simplifies  to  univariate  polynomial  splines.  For  a recent  review  of 
radial  basis  function  approximations,  see  [5]. 

This  note  is  about  applications  for  radial  basis  functions  and  other  approximation 
schemes  such  as  Sobolev-orthogonal  polynomials  and  more  general  Sobolev-orthogonal 
functions  to  the  numerical  solution  of  partial  differential  equations.  The  basic  ideas  stem 
from  the  theory  of  Sobolev-orthogonal  polynomials  ([13]),  and  in  this  paper  there  is  a 
remarkable  connection  developed  between  applications  of  Sobolev-orthogonality  with 
radial  basis  functions  (e.g.  [5]),  and  wavelets  are  mentioned  as  well  (e.g.  [8,  9]).  Sobolev- 
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orthogonal  polynomials  are  a device  to  extend  the  standard  theory  of  orthogonal  polyno- 
mials (see,  for  instance,  [12])  by  requiring  orthogonality  with  respect  to  non-selfadjoint 
inner  products  of  the  form 


(f,9) A=  [ f(x)g(x)dx  + 


for  a positive  parameter  A and  a suitable  interval  (a, 6),  a, b £ RU  (±oo).  The  (\x  in 
the  two  integrals  is  often  replaced  by  more  general  Borel  measures,  d^>,  say.  The  scheme 
which  we  want  to  discuss  in  this  short  article  is  one  of  spectral  type:  in  lieu  of  e.g. 
finite  element  spaces  as  underlying  piecewise  polynomial  approximation  spaces  for  the 
solution,  we  take  purpose-build  approximations  which  make  the  linear  systems  which  we 
need  to  solve  particularly  simple,  sometimes  even  diagonal. 

Therefore,  in  the  first  instance,  we  develop  a theory  of  applying  Sobolev-orthogonal 
polynomial  basis  functions  for  the  numerical  solution  of  partial  differential  equations  via 
a spectral  method.  Then  we  extend  this  idea  to  general  classes  of  radial  basis  function- 
type  methods,  where  shift-invariant  approximation  spaces  are  generated  with  Sobolev- 
orthogonal  basis  functions.  Due  to  the  introductory  character  of  this  paper,  our  dis- 
cussion is  restricted  to  relatively  simple  cases.  Our  presentation  is  illustrated  with  the 
one-dimensional  reaction-diffusion  partial  differential  equation. 

This  is  the  place  to  note  that  radial  basis  functions  have  found  a number  of  other 
applications  in  the  discretisation  of  PDEs.  Thus,  for  example,  Driscoll  and  Fornberg  [10] 
have  used  fast-converging  ‘flat’  multiquadrics  in  pseudospectral  methods,  while  Frank 
and  Reich  [11]  applied  radial  basis  functions  with  particle  methods  in  order  to  conserve 
enstrophy  in  the  solution  of  certain  shallow-water  equations.  Our  application  is  of  an 
altogether  different  nature. 


1.1  Examples  of  PDEs  and  Sobolev-orthogonality 

Consider  the  partial  differential  equation 

~~  = V (oVu)  + bu  + c,  (1.1) 


where  u — u(x,  t)  is  of  sufficient  smoothness  with  respect  to  x and  t,  x is  given  in  a cube 
V C Hd  (more  generally,  in  a finite  domain),  t > 0,  a = a(x)  > 0,  b = b(x)  and  c = c(x). 
We  impose  zero  Dirichlet  boundary  conditions.  The  stipulation  of  cube  as  a domain  and 
zero  Dirichlet  conditions  is  unduly  restrictive,  but  it  will  suffice  for  the  short  presentation 
in  this  paper  and  adequately  illustrate  the  main  novel  concepts  in  our  presentation.  In 
the  next  section,  we  shall  also  introduce  a nonlinearity  into  the  underlying  PDE. 

We  wish  to  approximate  the  solution  u(x,t ) as  a finite  linear  combination  of  the 
generic  form 

m 

u(x,t)  = ^ai(x)wi(t), 
i=i 

where  t is  normegatiye  and  x resides  in  the  domain.  In  the  sequel  we  shall  also  use 
expansions  into  infinite  series  with  l £ ZZ.  Thus,  a Galerkin  ansatz  (in  the  usual  L2  inner 
product  on  ]Rd  which  we  denote  by  (•,  ■)  in  contrast  to  the  specialised  Sobolev-inner 
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product  (•,  -)\  above)  gives 

m m m 

y\auak)w'l  =^y2(V(aVai),ak)wi  + Y ](bai,ak)wi  + ( c,ak ),  k = 1,2, . . . ,m. 

1=1  1=1  i=i 

Integration  by  parts  in  the  second  term  above  and  substitution  of  the  requisite  zero 
boundary  conditions  yield  the  alternative  formulation 


^(a«,  ak)w[  = - (aVa;,  Vak)  wt  + ]P(6a;,  a*)^  + (c,  ak ), 


1, 2, . . . , m. 


We  solve  the  ODE  system  (1.2)  with  respect  to  t,  for  example  with  the  backward  Euler 
scheme  (we  use  backward  Euler  for  the  sake  of  simplicity,  but  it  should  be  noted  that 
the  same  analysis  applies  to  any  implicit  multistep  method,  because  our  use  of  Sobolev- 
orthogonality  is  only  linked  to  the  implicitness  of  the  solution  method) 

ly?*1  = wf  + AtFi(wn+1),  n £ 2Z+,  J = l,2,...,m,  (1.3) 

where  the  function  f)  is  given  implicitly  by  the  equations  (1.2)  and  where  w"+1  in 
the  expression  above  is  the  vector  with  components  u>"+1,  l — 1, 2 ,m.  Let  us  now 
multiply  expression  (1.3)  by  ( ai,ak ) and  sum  up  for  l = 1,2, . . . ,m.  Then,  exploiting 
(1.2),  a little  algebra  yields 

m [1  - At6(x)]  Q/(x)aA-(x)dx  + At  J a(x)VTQi(x)Va;.(x)dx|  w"+1 

= Y]  I ai(x)afc(x)dxu;/n  + / c(x)afc(x)dx.  (1.4) 

JZ (Iv  Jv 

The  connection  with  Sobolev-inner  products  is  clear.  Indeed,  let  us  now  choose  the  set 
Wm,n  ■:=  {wi,W2,  ■ ■ . , wm } as  a set  of  functions  that  are  orthogonal  with  respect  to  the 

O 

homogeneous  Sobolev  Hd, 2 inner  product  (see,  e.g.,  [13]) 

(/)  9) At  ■■=  [ [1  - At&(x)]/(x)g(x)dx  + At  f o(x)VT/(x)V(/(x)dx  (1.5) 
Jv  Jv 

(this  of  course  requires  that  At6(x)  < 1,  hence  may  restrict  in  a minor  way  the  choice  of 
the  time  step  At).  Further  below  we  shall  also  use  infinite  sets  W instead  of  the  finite 
set  W„,i7l.  It  is  important  to  note  that  in  general  the  Sobolev  inner-product  depends 
upon  the  step  size.  Subject  to  this  formulation,  the  linear  system  (1.4)  diagonalises  and 
its  numerical  solution  becomes  trivial.  We  turn  now  to  a more  elaborate  example  in  the 
next  subsection,  namely  the  reaction-diffusion  equation. 

1.2  Reaction-diffusion  as  a paradigm  for  nonlinear  PDEs 

Let  us  consider  the  nonlinear  partial  differential  equation 

^ = V(aV«) + /(«),  (1.6) 
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where  otherwise  all  the  quantities  are  as  in  (1.1),  including  the  boundary  conditions. 
Suppose  that  an  approximation  un  to  tt(x,  nAt)  is  available  at  all  the  spatial  grid  points. 
We  commence  by  interpolating  un  to  requisite  precision  by  some  function  v.  Thus,  v is 
defined  throughout  the  cube  V and  coincides  with  un  at  the  grid  points.  This  allows  us 
to  linearise  the  source  function  / about  un,  the  outcome  being 

Ou 

= V(aVu)  + c + bu  + g(u),  (1.7) 

where 


b(x)  = />(x)), 

c(x)  = f(v(x.))  - f'(v(x))v(x), 

£/(x,u)  = f(u)  - f(v(x.))  - f(v(x))[u-v(x)]. 

Note  that 

g(x,u)  = 0(\u-v\2). 

We  can  now  solve  the  nonlinear  system  (1.7)  by  functional  iteration,  i.e.  by  letting  as  a 
start 

w^+1’0  = wf,  l = 1,2, ...  ,m, 


and  recurring,  employing  the  inner  product  (1.5), 


J2(ai,ak)Atw?+1'j+1  (1.8) 

i=i 

^(cq,  ak)w f + (g  ( ■ , , ak ) , k = 1, 2, . . . , m, 


i=i 


i=i 


for  j € 2Z+. 


If,  as  in  the  previous  subsection,  we  choose  Wm  so  as  to  diagonalise  the  linear  sys- 
tem, each  step  of  (1.8)  becomes  relatively  cheap.  Hence  this  approach  might  offer  a 
realistic  means  to  derive  spectral  approximation  to  nonlinear  PDEs.  Indeed,  a special 
one-dimensional  case  can  be  treated  straightforwardly  and  it  is  presented  in  the  sequel. 


1.3  The  one-dimensional  case  using  polynomial  splines 


Let  (1.1)  be  given  in  one  space  dimension  and  without  source  terms,  whence  it  becomes 
the  familiar  diffusion  equation  with  variable  diffusion  coefficient, 


du  _ d f du\ 
dt  dx  \ dx ) 


Thus,  provided  that  0 < x < 1 and  t nonnegative,  we  require  the  ‘usual’  Sobolev 
orthogonality  [13]  with  respect  to  the  inner  product 


(/)  9)  At  = (f,g)  1 


f 


f(x)g(x)d<p{x)  + 


f f'{x)g'{x) dif}(x), 
Jo 
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= 1 - A tb, 


where 

iM  = 1-Atb,  Mh- Ma. 
dx  dx 

We  emphasise  again  the  dependence  of  the  Sobolev-inner  product  on  the  step  size.  Taking 
the  approach  of  the  previous  subsection  as  our  point  of  departure,  an  obvious  option  is 
to  use  Sobolev-orthogonal  polynomials.  An  alternative  approach  which  can  be  worked 
out  explicitly  and  which  we  wish  to  demonstrate  in  this  subsection,  is  to  use  univariate 
polynomial  spline  approximations.  It  has  the  advantage  of  being  more  amenable  to  a 
generalisation  to  several  space  dimensions. 

We  suppose  that  the  unit-interval  [0, 1]  is  divided  into  N intervals  of  length  h :=  jj 
and  consider  a piecewise-quadratic  basis  of  continuous  functions  «i,  «2i  • • • , Sn  such  that 

^[x  - (l  - l)/i]  + ai(x  - lh)[x  - (l  — l)/i],  (Z  — 1 )h  < x <lh, 

Si(x ) :=  < j^[(l  + l)/i  — x]  + fa(x  - lh)[x  — (l  + 1 )h],  Ih  < x < (l  + 1 )h, 

k 0,  |x  - lh\  > h. 

Clearly,  si  is  a continuous,  C[0, 1]  cardinal  function  of  Lagrange  interpolation  at  the 
knots  (hence,  a quadratic  spline  with  double  knots,  cf.,  Powell  [16],  the  added  degree  of 
freedom  taken  up  by  the  requirement  of  Sobolev-orthogonality).  Next,  we  need  just  to 
impose  Sobolev  orthogonality,  and  solve  for  the  coefficients  eq  and  0i . This  is  equivalent 
to  the  requirement  that 

(si,si+i)&t  = 0,  l = 1,2, . . . ,N  - 1. 

In  the  special  case  a(x)  = 1,  6(x),c(x)  = 0,  we  have  ip(x)  = x,  ip(x)  = A tx  and 


(s;,s*+l)At 


+ At  J ^ — + 2ai+ix  — ai+ih'j 


fax(x  — h)  dx 


- + 2 ai+lx  - ai+ihj  ( ■ --  + 2 fax  - fahj  dx 


= h£[t  + - lk]  • [l  - e - Pih2Z(  1 - 0]df 

At 

— y (!  + 2 ai+1£  - am)(l  + fa-  2/?,£)d£ 

, [71  h2 , h 4 \ At  / 1 

" h [\6  ~ 12{ai+l  + Pl)  + 30 ai+lf3lJ  + V ( 1 + 3 


1 + oQ*  + lA 

& J 


Let  n = At/ h2  be  the  Courant  number.  Since  we  have  two  degrees  of  freedom  for  each 
l and  because  each  equation  is  otherweise  independent  of  l,  we  may  fix  q:  = a;  s fa. 
Then,  letting  a :=  h2a,  requiring  (s;,S(+i)a*  — 0 is  equivalent  to 

5 - 5d  + &2  + 10/ic*2  — 30/r  = 0 (1.9) 


(10  n + h4)a2  — 5 h2a  + 5 — 30  p = 0. 
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We  wish  to  solve  this  quadratic  equation  for  a for  a suitable  range  of  Courant  numbers. 
Indeed,  the  equation  (1.9)  has  two  real  solutions  a for  every  p > | if  h is  small  enough, 
since  its  discriminant  is 

(120/z  + 5 )hA  + 1200 /i2  - 200 p. 

In  the  case  p = | each  s;  reduces,  upon  the  choice  of  a = 0,  to  a chapeau  function. 
Otherwise  we  obtain  a = 0(1).  We  may  give  up  a small  support,  characteristic  of 
spline  functions  (which,  anyway,  is  of  marginal  importance,  since  we  do  not  solve  linear 
systems!).  This  is  a case  discussed  in  the  next  section.  Another  obvious  alternative  is  to 
construct  an  orthogonal  basis  from  chapeau  functions.  This,  however,  is  easily  seen  to 
be  identical  to  the  LU  factorization  of  the  standard  FEM  matrix 

0 0 0 •••  0 ' 

| | 0 ° •••  ° 

0 | | § 0 •••  0 

° •••  0 i | I 0 

0 ...  0 0 I I i 

U ...  0 0 § 1 1J 


2 Applications  of  radial  basis  functions  and  wavelets 


2.1  Sobolev-orthogonal  translates  of  a radial  basis  function 

In  this  section,  we  wish  to  develop  a more  general  approach  employing  the  concepts 
of  wavelets  and  radial  basis  functions  and  employ  shift-invariant  spaces  of  approxima- 
tions for  our  spectral  methods.  We  begin  by  giving  up  the  compactness  of  the  domain 
V and  work  on  the  entire  real  line  instead.  For  this,  we  shall  demonstrate  the  use  of 
Sobolev-inner  products  and  shift-invariant  spaces  and  concentrate  solely  on  this  part 
of  the  analysis  in  the  present  article.  So,  in  particular,  the  set  W above  is  of  the  form 
{</->(■  — nh)  | n e 2Z}.  In  the  sequel  we  shall  add  several  remarks  about  how  to  find 
compactly-supported  <f>  that  allow  the  treatment  of  partial  differential  equations  on  com- 
pact domains.  We  remark  that  n is  no  longer  used  for  the  time-steps  in  the  differential 
equation  solver  but  for  the  shifts  of  the  radial  functions. 

To  start  with,  we  wish  to  find  a function  <j>  € H2(IR),  where  H2(1R)  is  a non- 
homogeneous  Sobolev  space,  such  that  for  a positive  constant  A and  positive  spacing 
h it  is  true  that 

/OO  pOO 

<j>{x)<p(x  - hn)  dx  + X <f>' (x)(j>' (x  - hn)  dx  = <50n)  n&7L.  (2.1) 

-oo  J — OO 


We  multiply  both  left-  and  right-hand-side  of  the  general  pattern  (2.1)  by  exp(i0n)  and 
sum  over  n € ZZ, 


°°  ( r oo  poo 

exp(i  On)  < / <j>(  x)<f{x  - hn)  dx  + X <j>’ (x)</>' (x  - hn)  dx  > = 1,  9 € [-7T,  ir\. 

%=- oo  U-oo  J- oo  J 

(2.2) 
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In  order  to  be  able  to  exchange  summation  and  integration  and  apply  the  Poisson  sum- 
mation formula  (Stein  and  Weiss  [17],  p.  252)  we  make  a number  of  assumptions.  The 
version  of  the  Poisson  summation  formula  that  we  wish  to  use  states  that  for  a univariate 
function  / with 

1/0*01 = °((i + Mr1-') 

and 

\f(x)\=o{(i+\x\r1-^ 

and  positive  e,  the  following  equality  holds  (note  that  the  first  bound  in  the  above  implies 
existence  and  continuity  of  the  one-dimensional  Fourier  transform) 

OO  OO 

f(9  + 2itj)  = ]T  exp(i 6j)f{j). 

jf=  — OO  j~  — OO 

Specifically,  we  assume  that  the  following  three  decay  estimates  hold: 

\<j){x)\  < c(  1 + |t|)_1_£, 

10' 0*01  < c(l  + |t|)_1_£, 


' and 

10(01  < c(l  -P  |^|)-3-e, 

where  c is  a generic  positive  constant,  e > 0,  4>  denotes  the  Fourier  transform  and  we 
demand  the  faster  rate  of  decay  in  the  last  display  because  we  shall  later  require  summab- 
ility  of  translates  of  the  Fourier  transform  multiplied  by  the  square  of  its  argument.  Note 
in  particular  that  the  first  decay  condition  renders  the  Fourier  transform  continuous 
and  well  defined. 

An  example  for  a function  <j)  that  satisfies  the  three  decay  conditions  above  is  the 
second  divided  difference  of  the  multiquadric  radial  basis  function  [4]  vV2  + C2  that  is 

<j>(x)  = ly/(x-  l)2  + C2-  \/x2  + C2  + ±y/(x  + l)2  + C2. 

Here,  C is  a positive  constant  parameter.  The  above  function  decays  cubically  [4]  and 
its  Fourier  transform  even  decays  exponentially  due  to  the  exponential  decay  of  the 
modified  Bessel  function  K\  [1]  that  features  in  the  generalised  Fourier  transform  of  the 
multiquadric,  here  stated  only  in  the  one-dimensional  case, 


(cf.  Jones  [14]). 

Once  summation  and  integration  are  interchanged,  (2.2)  becomes 


/OO  00 

4>{x)  E exp(i9n)(j)(x  - hn)  da; 

•OO  „ — 
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/OO  00 

4>'{x)  E exp(i 6n)(j)'{x  — hn ) dx  = 1,  9 £ 

-OO  


[— 7T,  7r],  (2.3) 


or,  applying  the  Poisson  Summation  Formula  (Stein  and  Weiss,  [17],  p.  252) 

/OO  ’ oo  . 

<p[x)  ^ exp^i/i-1a;(0  + 2'irn)\(f>(h~1(6  + 27m)  j da;  + iA/i-1  (2.4) 

■°°  n=— 00 

/OO  00 

^'(a:)  ^ expfih_1a;(0  + 27rn)'j  (9  + 2irn)4>(h~1(0  + 27rn))  da:  = h, 

where  9 € [— 7r,  it\.  Because  (p  vanishes  at  infinity,  integration  by  parts  of  the  second  term 
of  (2.4)  gives 

/oo  00 

<j)(x)  ^ exp(^ih~1x(6  + 27rn)^$(ji~1(6  + 2nn)Sj  dx 

■°°  n=— 00 

^ /»oo  00 

+ V2  / 0(x)  exp(i/i_1a;(0  + 2itn)){8  + 2nn)2(p(h~l(9  + 27rn))da; 

-'-oo  „=_oo 

00 

= y 4>(h~l{6  + 2rn))d)(^—h~l  (9  + 27rn)^  l + \h~2 (6  + 2nn)2  = h. 

n=— 00 

Since  (p  is  real,  £)  = </>(£),  and  this  implies 

OO 

[^(/i-1^  + 27m))[2(l  + A h~2(0  + 27m)2)  = h,  0e[-7r,7r].  (2.5) 

n~~  00 

This  is  our  condition  that  leads  to  the  required  Sobolev-orthogonality.  In  summary,  we 
have  established  the  following  theorem. 

Theorem  2.1  If  the  decay  conditions  on  <p,  as  stated  above,  hold  in  tandem  with  the 
expression  (2.5),  then  the  required  orthogonality  condition  (2.1)  is  satisfied. 

We  note  that,  if  we  are  given  a ip  such  that 


{6  + 27rn)^ 

1 = h,  9 € [— 7T,  7r], 

(2.6) 

m :=  - 

A 

HO 

A + A£2 

(2.7) 

satisfies  (2.5).  This  expression  can  be  used  to  derive  an  explicit  transformation  which 
takes  a ip  that  satisfies  (2.6),  into  a (p  satisfying  (2.5),  although  its  practical  computation 
may  be  nontrivial.  Indeed,  by  the  Parseval-Plancherel  theorem  [17],  we  get  the  useful 
identity 


^(X)  = w^/  ^(x~y')K°  (^)  d?/> 


206  M.  Buhmann,  A.  Iserles,  and  S.  P.  N0rsett 

which  is  a convolution  and  whose  Fourier  transform  is  therefore  (2.7)  (cf.,  for  instance, 
Jones  [14]).  In  (2.8),  Kq  is  the  0th  modified  Bessel  function  (Abramowitz  and  Ste- 
gun  [1])  which  is  positive  on  positive  reals  and  satisfies  K0(t)  ~ — log  t near  zero  and 
K0(t)  ~ y/7r/(2f)e"t  for  large  t,  similar  to  the  asymptotics  We  have  used  before  for  the 
K\  modified  Bessel  function.  Hence,  by  a lemma  in  [7],  see  also  (Light  and  Cheney  [15]) 
<!>  decays  algebraically  of  a certain  order  if  y)  does.  Moreover,  because  1 / %/ 1 + Ax2  is 
positive,  integer  translates  of  <p  are  dense  in  L2,  say,  provided  that  this  is  the  case  with 
, integer  translates  of  ip  [18]. 

In  some  trivial  cases  we  may  evaluate  the  integral  (2.8)  explicitly,  for  instance  for 
ip(x)  — cosx,  where  the  integral  is  again  a constant  multiple  of  the  cosine  function 
(Abramowitz  and  Stegun  [1]).  Otherwise,  the  smoothness  and  fast  exponential  decay  of 
the  modified  Bessel  function  can  be  used  together  with  a quadrature  formula. 

We  may  now  use  the  translates  of  such  Sobolev-orthogonal  functions  in  the  spectral 
approximation  of  a PDE  as  above,  letting  W :=  {<f>(-  — nh ) | n € 2Z }. 

An  example  of  a function  ip  that  satisfies  (2.5)  is  simply  the  characteristic  function 
scaled  by  h of  the  interval  [— hn,  hit].  In  that,  case,  \ip(x)\  decays  like  l/|x|.  In  fact,  any 
ip  that  satisfies  |t/>(£)|  < c(l  + |£|)_1/2_e  for  positive  e can  be  made  to  satisfy  (2.6)  by 
subjecting  it  to  the  transformation 


see  for  instance  (Battle  [2]).  If  ip  is  compactly  supported  then  the  transformed  y'j  will 
not  necessarily  be  compact  supported  but  decay  exponentially  [6] . 

In  order  to  find  a class  of  examples  of  compactly  supported  ip  that  satisfy  (2.6),  see 
Daubechies  [8]  for  her  compactly  supported  scaling  functions  ip  which  are  fundamental 
for  the  construction  of  Daubechies  wavelets.  For  example,  the  following  conditions  are 
sufficient  for  ip  which  shall  be  defined  by  its  Fourier  transform  to  satisfy  (2.6)  for  h — 1 
(other  h can  be  used  by  scaling): 

OO  £ 

'to'II'W 

i=l 

where,  for  some  suitable  coefficients  hk, 

2N-1 

MO  = £ hke 

fc= o 

has  to  satisfy  h{ 0)  = 1,  h{ir)  = 0,  and 

IM0l2+1M£  + 7r)|2  = 1.  £e[-7T,7r]. 

For  the  construction  of  such  h,  see  [8] . Compactly  supported  basis  functions  are  import- 
ant to  approximate  the  numerical  solution  of  a PDE  as  in  the  above  example  defined  on 
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a compact  V.  Moreover,  any  if  with  the  aforementioned  decay  property  can  be  made  to 
satisfy  (2.5)  by  the  transformation 


4>(Q  ” 


y/hip{  Q 

oo 


'Y,  | ip(f  + h 127rn)|2(l  + A(£  + /i_127m)2 


(2.10) 


They  can  also  be  found  by  applying  the  transformation  (2.10)  and  using  the  transform- 
ation (2.9)  as  well. 

We  note  finally,  that  for  instance,  when  ip  is  a B-spline  then  its  translates  are  dense 
in  L2  if  we  allow  h to  become  arbitrarily  small  (see,  for  instance,  Powell  [16])  and  the 
last  section  of  this  paper).  .' 

2.2  Sobolev-orthogonal  translates  of  a function  in  higher  dimensions 

Applying  the  approach  of  the  previous  subsection  to  the  Sobolev  inner  product 

/ /(x)ff(x)  dx  + A f VT/(x)V#(x)dx, 

J Rd  J Rd 


the  outcome  is  the  orthogonality  condition 


\H^1(0  + 2irn))\2(l  + Xh-2\\e  + 2nnf)  = hd,  0e[-7r,7r]d,  (2.11) 


which  replaces  (2.5).  We  are  now  also  interested  in  the  more  general  case  of  Sobolev-type 
inner  products 


/(x)s(x)/i(x)  dx  + A f VT/(x)V^(x)z/(x)  dx, 

J RA 


where  the  weights  /x  and  v are  positive.  Here  the  orthogonality  condition  becomes  more 
complicated.  Specifically,  it  is 


Y l{0  + 1(0  + 27rn)^ 

n e7Ld 

+ \h~2<pv(h~l(9  -)-  27rn)^„^h_1(0  + 27rn)j  =hd,  6e[— n,Tr]d, 

where 

0/t  ;==  4>  * \/m> 

4>v  :=  (INI  x $)  * y/u, 

and  * denotes  continuous  convolution,  used  as  in  (2.8),  where  ip  is  convolved  with  a 
modified  Bessel  function. 

2.3  Error  estimates 


We  can  offer  error  estimates  for  the  Sobolev-orthogonal  bases,  firstly,  in  the  case  when 
(p  is  a univariate  spline  of  fixed  degree  m,  say,  with  knots  on  h7L , and,  secondly,  in  the 


M.  Buhmann,  A.  Iserles,  and  S.  P.  N0rsett 


208 

case  when  0 is  a linear  combination  of  translates  of  the  radial  Gauss  kernel 

e-"2*2/2,  x e P, 

along  h7L.  In  the  former  case  it  is  known  that  the  uniform  approximation  error  to  a 
sufficiently  smooth  function  from  the  linear  space  spanned  by  </>(•  - nh),  n £ 77,  is  at 
most  a constant  multiple  of  hm+1  ([16]).  We  have  already  mentioned  that  we  require 
A — 0(h2),  therefore  it  can  be  deduced  by  twofold  integration  by  parts  that  the  Sobolev 
error  is  indeed  0(hm+1).  This  can  be  generalized  in  a straightforward  way  to  higher 
dimensions  by  tensor-product  B-splines. 

Our  L2(1R)  error  estimates  can  be  carried  out  as  follows:  Let  / be  a band-limited 
function,  that  is,  one  with  a compactly-supported  Fourier  transform,  which  satisfies 
such  assumptions  that  imply  that  the  best  least-squares  approximation  using  a Sobolev 
inner  product 

OO 

Sh(x)  = (f,<p(-  - hn))yh<f>(x  - nh),  x £ ]R,,  (2.12) 

n=— oo 

is  well  defined.  For  instance,  we  may  require  that  ( < oo,  as  well  as  sufficient 
decay  of  the  radial  basis  function  <f>,  i.e. 

\4>(r)\  < c(l  -f-  |w|)-1-e, 

\(t>\r)\  < c(l  + M)-1~£, 

l^(r)l  < c(l  -F  Irl)-1-^ 

for  a positive  e.  Here  (-,-)x,h  is  the  Sobolev  inner  product  which  we  study  in  this  note 
and  it  is  helpful  to  emphasise  its  dependence  on  h in  the  subscript.  We  begin  with  the 
piecewise  polynomial,  i.e.  spline,  case.  Hence,  let  <j>  be  from  the  space  of  splines  of  degree 
m with  knots  on  h7L  such  that  its  translates  are  Sobolev  orthogonal. 

Theorem  2.2  Subject  to  the  assumptions  of  the  last  paragraph,  we  have  the  error 
estimate 

\\sh-f\\2  = 0(hm+l),  h-+  0.  (2.13) 

Proof:  We  shall  establish  in  the  course  of  this  proof  an  error  estimate  for  the  first 
derivative  of  the  error  function  in  (2.13),  so  that  an  order  of  convergence  can  also  be 
concluded  for  the  norm  associated  with  our  Sobolev  inner  product.  Indeed,  because  the 
Fourier  transform  is  an  L2(R)  isometry,  we  may  prove  (2.13)  by  considering 

Wh  ~ fh  (2-14) 

instead  of  the  left-hand  side  of  (2.13).  The  Fourier  transform  of  (2.12)  is 

OO 

h{0)  = ]T  (/,  0(.  - nh))x  he~ienhm,  e e H. 

n=— oo 

The  absolute  convergence  of  the  above  is  guaranteed  by  the  decay  conditions  on  <j>. 
Hence  the  square  of  (2.14)  is,  by  the  Parseval-Plancherel  Formula  and  periodisation  of 
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de 


the  integrand  with  respect  to  9, 

/OO  00 

f(0)  - E (f,<P(--nh))Xhe-lShn$(9) 

‘°°  n=—OC 

/OO  I OO  «00 

/»-  E / /(0^)e^n(l  + A^2)d^e 

-ool  n=-ooJ-°° 

/ir/h  °°^ 

'E  f {9  + 2nk / h)  — (j)(6 -\- 2itk / h) 

"*/hk=-oc 


-Whn 


m 


dd 


OO  -OO 

E / hm0^nh(l  + A£2)  d^e-10™'1 

_ ^ J — 00 


d9. 


(2.15) 


The  (1  + A£2)  term  in  the  above  comes  from  the  derivative  in  the  Sobolev  inner  product 
and  Fourier  transform.  Because  / is  band-limited,  for  small  enough  h (2.15)  assumes  the 
form 


d9. 


/Tf/h  00  00  poo 

E f(e)60k-4>(e  + 2nk/h)  E / mm)^nh( l + \e)dte-ienh 

-7r/,1fc=- 00  n=-ooJ-°° 

(2.16) 

Using  again  the  band  limitedness  of  /,  together  with  the  Poisson  Summation  Formula, 
(2.16)  can  be  brought  into  the  form 

/7T /h  °°  I 

E \mSok-f>(0  + 2nk/h) 

■ir/h.  1 

|2 


'lh  k=- 00 

OO 

h 


1 ^ , 

X-  E /(#  + 2im/h)(t>(0  + 27rn/h)  ^1  + A(0  + 27rn//i)2J 


d(9 


/ir/h  00 

E |/(*)*0k  - + 27rfc/h)/(0)<£(0)(l  + A6>2)|  de. 

■*/hk=- OO 


(2.17) 


In  the  case  when  (f>  is  in  the  aforementioned  spline  space,  it  can  be  expressed  as  the 
inverse  Fourier  transform  of 


m 


Vhr{£) 


\Jl2n=-oo  |r(£  + ft_127m)|2(l  + A(£  + h-127m)2) 


(6R, 


(2.18) 


where  f(f  ) = £_m_1.  This  follows  from  (2.5)  and  from  the  fact  that  all  splines  from  our 
space  are  linear  combinations  of  integer  translates  of  r(x)  :=  \x\m,  whose  generalised 
Fourier  transform  is  a multiple  of  £~m~1  [14].  Since  any  constant  factors  in  front  of 
the  function  £~m~1  in  f cancel  in  the  expression  for  <j>  above,  we  have  ignored  them 
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straightaway.  Substituting  (2.18)  into  (2.17),  we  get  the  integral  over  [-w/h,  n/h]  of 


E fm o*--55 f(e+hr.2nkM9) im  + \e2)  . (2.19) 

'=-°°  E \r(0  + h-127m)\2{l  + Xid  + h-^rcn)2) 


Considering  (2.19)  for  each  m separately,  it  follows  from  (2.19)  and  from  f(£)  = £-m-1 
that  our  claim  is  true.  Indeed  for  the  sum  over  all  terms  with  k ^ 0,  it  is  evident  that 
we  obtain  a factor  of  h2m+2  from  the  numerator,  because  the  denominator  is  periodic, 
containing  one  term  independent  of  h,  and  the  nonvanishing  expression  h~12nk  in  the 
argument  of  f(8  + h~12irk)  guarantees  f{8  + h~12nk)  ~ hm+1  due  to  r(£)  = Of 

course,  the  squares  then  taken  provide  the  h2rn+2  instead  of  hm+1. 

On  the  other  hand,  for  k — 0,  we  have  for  small  enough  h 

f(a) if(g)i2(i  + Ae2)/» 2 

S^=-oc  1^(0  + h_127rn)|2(l  + \(8  + h~12rrn)2) 

^njiol^  + h-^Tm^il  + Xie  + h-^Tm)2)  2 

1 + + h_127rn)|2(l  + X(8  + h-12irn)2) 

which  is  also  0(h2m+2),  as  required,  because  the  numerator  provides  an  0(h2m),  ac- 
cording to  the  rate  of  the  decay  of  r and  the  power  of  h in  its  argument.  This  is  then 
squared  to  provide  0(h4m)  = 0(h2rn+2). 

As  for  the  derivatives,  one  only  has  to  multiply  the  Fourier  transform  of  the  error 
function  in  (2.14)  with  8,  and  we  get  the  same  error  estimate  by  multiplying  the  integ- 
rands in  all  the  following  integrals  with  \8\2.  □ 

The  same  analysis  remains  valid  when  considering  integer  translates  of  the  Gauss 
kernel  e-7  x /2  in  order  to  form  (p.  In  this  case  we  make  use  of  the  fact  that  the  Gauss 
kernel  has  a Fourier  transform  which  is  a multiple  of  e~x  A2a  ).  \ye  put  this  instead  of 
f into  (2.19),  and  we  then  get  arbitrarily-high  orders  of  convergence  from  (2.14)  as  long 
as  we  take  7 = 0(h),  see  also  [3].  For  this  choice  (p  is  exponentially  decaying,  whereas 
for  splines  of  degree  m we  merely  get  algebraic  decay  at  infinity  of  order  -m  - 1. 
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